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Coupling light to Rydberg states of atoms under conditions of electromagnetically induced transparency (EIT) 
leads to the formation of strongly interacting quasi-particles, termed Rydberg polaritons. We derive a one¬ 
dimensional model describing the time evolution of these polaritons under paraxial propagation conditions, 
which we verify by numerical two-excitation simulations. We determine conditions allowing for a description 
by an effective Hamiltonian of a single-species polariton, and calculate ground-state correlations by use of the 
density matrix renormalization group (DMRG). Under typical stationary slow-light EIT conditions it is difficult 
to reach the strongly interacting regime where the interaction energy dominates the kinetic energy. We show that 
by employing time dependence of the control field the regime of strong interactions can be reached where the 
polaritons attain quasi crystalline order. We analyze the dynamics and resulting correlations for a translational 
invariant system in terms of a time-dependent Luttinger liquid theory and exact few-particle simulations and 
address the effects of nonadiabatic corrections and initial excitations. 

PACS numbers: 32.80.Ee,32.80.Qk,42.50.Gy 


I. INTRODUCTION 


Recently there has been growing interest in using Rydberg 
gases |T1 to tailor strong and nonlocal nonlinearities for pho¬ 
tons. The strong interactions between Rydberg states in a gas 
of ultra-cold atoms can mediate strong nonlocal interactions 
for light fields propagating under conditions of electromagnet¬ 
ically induced transparency (EIT) in such a medium |2}{T(j). 
For instance, the formation of a small blockaded volume lead¬ 
ing to antibunching of photons as well as pronounced bunch¬ 
ing of photons in a Rydberg gas have been demonstrated ex¬ 
perimentally HUE). Moreover, first steps have been taken 
towards building single photon logic devices, such as single 
photon switches EHE). 

Rydberg polaritons exhibit potentially strong and nonlo¬ 
cal mutual interactions making them interesting candidates 
to study many-body effects in the regime of strong correla¬ 
tions, where the interaction energy dominates the kinetic en¬ 
ergy. In the present paper we consider the propagation of light 
in Rydberg EIT media under paraxial conditions in a one¬ 
dimensional setting, as can be realized experimentally by, e.g., 
cigar-shaped atomic ensembles or sending photons through 
hollow-core optical fibers filled with Rydberg atoms HUE). 
In a previous letter |fT9l we argued that the propagation of 
photons inside Rydberg gases under EIT conditions can be 
described in terms of interacting quasi-particles, named Ryd¬ 
berg polaritons. We showed that a regime of strong correla¬ 
tions, where the ratio between interaction energy and kinetic 
energy becomes large, can be reached by dynamically turning 
photons into stationary Rydberg excitations ll20l . 

In the present paper we extend these studies by deriving an 
effective dark-state polariton model including leading-order 
corrections valid for sufficiently large separation of Rydberg 
polaritons. This regime allows for a perturbative treatment of 
the coupling between bright- and dark-state polaritons. We 
derive a Lindblad master equation for the Rydberg dark-state 
polaritons where the bright-state polaritons act as a Marko¬ 
vian reservoir. We verify the model and address the prepara- 
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Figure 1. a) Atomic level scheme with input probe field £ and control 
field fl. b) Sketch of a possible experimental realization: one dimen¬ 
sional setup with copropagating photons focused inside an elongated 
atomic cloud of Rydberg atoms and control field Q. 


tion of an initial state by employing numerical two-excitation 
wave-function simulations. We derive conditions for when 
the Master equation can be reduced to the effective Hamil¬ 
tonian introduced in lfl9) . Furthermore, we show that for an 
excitation density smaller than one per blockade volume and 
a transversal beam diameter less than the blockade radius the 
paraxial propagation can be described by a one-dimensional 
model. 

The paper is organized as follows. In Section II we intro¬ 
duce the microscopic model of photons coupled to interacting 
three-level atoms. We map the model to the polaritonic basis 
and derive a one-dimensional Born-Markovian master equa¬ 
tion for the dark-state polaritons by treating the bright-state 
polaritons as a reservoir. Calculating interaction-mediated 
scattering amplitudes between different transverse modes we 
derive conditions for an effective one-dimensional treatment 
of the problem. We furthermore analyze conditions when 
the master equation reduces to an effective Hamiltonian de¬ 
scribing a unitary time evolution. In section III we employ 
numerical two-excitation wave-function calculations to ver¬ 
ify the model and determine the initial state after sending 
a light pulse into a Rydberg gas. In section IV we calcu¬ 
late ground state correlations of the effective Hamiltonian by 
use of density-matrix renormalization group (DMRG) calcu¬ 
lations and a Luttinger liquid approach. We show that un- 
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der continuous driving conditions and large optical depth per 
blockade distance the regime of strong interactions cannot be 
entered since the Rydberg blockade prevents to reach the po- 
lariton densities necessary for the interaction energy to over¬ 
come kinetic contributions. However, in section V we show 
that it is possible to enter this regime by dynamically slow¬ 
ing down Rydberg polaritons. In particular we consider the 
storage of interacting photons inside a Rydberg medium. We 
apply numerical wave function simulations to show that the 
losses during storage remain finite and use a time-dependent 
Luttinger liquid approach to compute time-dependent density 
correlation functions. 


n. RYDBERG POLARITONS 

In this section, we discuss the paraxial propagation of weak 
light pulses consisting of few photons in an atomic medium 
under EIT conditions with Rydberg interactions. 


A. Light-matter coupling 

Let us consider a system of a weak quantized probe field 
E = yJ ^ E.£ e -i(u P t-k P z ) _|_ h a propagating through an en¬ 
semble of non-interacting three level atoms as sketched in 
Fig. B b) with the level structure in ladder configuration as 
in ED, and as drawn in Fig. E a). The operators £ , £' are 
slowly varying envelope operators, which obey bosonic com¬ 
mutation relations [£(r), £ 5(r')] = S(r — r') and uj p and k p 
denote carrier frequency and wave vector of the probe field, 
respectively. The atoms are composed of a ground state | < 7 ), an 
intermediate excited state |e) and a metastable Rydberg state 
| r). Here we neglect dipole-dipole interactions of Rydberg 
states for a moment, which we will reintroduce later. The field 
£ is coupled to the coherence between atomic states \g) and 
\e) with single-photon detuning A = w eg — w p . The coher¬ 
ence between states \e) and |r) is driven by a classical control 
field Q with carrier frequency uj c and detuning A c = w re — w c . 
We denote the resulting two-photon detuning by <5 = A + A c . 
The atomic state \e) is assumed to be subject to spontaneous 
decay with rate 7 . 

The atomic polarization is microscopically described by 
spin flip operators |/r) (v\ of individual atoms. By averaging 
these over a small volume centered at a position r and con¬ 
taining N r 1 atoms we define continuous, coarse-grained 
atomic spin flip operators a pl , at position r, 

1 Nr 

<V(r) = — (1) 

r j =1 

which fulfill the commutation relations [< 5 ‘ Q: / 3 (r), <5"^^(r')] = 
^S( r - r') [6pfj,a a „( r) - , where the density n of 

the atoms is assumed to be homogeneous. Transforming to a 
frame rotating with the atomic frequencies and performing the 
rotating wave approximation, the atom-light coupling Hamil¬ 


tonian of this system can be written as (h = 1 ) 

H = n J d 3 r jAcr ee (r) + <5a> r (r) j 

- jn<T re (r) + gy/ni(r)a eg (r) +h.a.j, (2) 

where g = d ge ^/u} ge /2He 0 denotes the coupling strength of 
the electric field £ to the atomic transition | g) — |e), with d ge 
being the atomic dipole matrix element. 

We derive Heisenberg-Langevin equations of motion for 
the slowly varying field £ and atomic operators taking into 
account the spontaneous decay rate 7 of the intermediate 
level l22l . Assuming the probe field to be weak compared to 
the control field we can treat the equations in linear response 
with respect to g £, leading to 

dtCTgr — 'l&CJtrY ~\~ 

® „ (3) 

d t i T ge = -Td-ge + igsjn£ + i£ld gr + F ge , 

where T = 7 + iA. F ge denotes a Langevin noise opera¬ 
tor Il22l . which is b-correlated in time and space with van¬ 
ishing expectation value that needs to be added to preserve 
the commutation relations. As the noise is related to the pop¬ 
ulation of the excited state which we set cr ee = 0 in linear 
response, the Langevin operators can be neglected. 

In the following we want to consider the case of two-photon 
resonant driving, i.e., <5 = 0. We note, however, that taking 
interactions between atoms excited to level |r) into account 
will lead to an effective space-dependent two photon detuning. 


B. Paraxial light propagation 

The dynamics of the probe field is described by the trun¬ 
cated paraxial wave equation, 

(I + 4 - 'W r 71 ) ^ (r ’ () = '9VSWM). (4) 

We assume a cylindrical symmetry of the setup and decom¬ 
pose the probe field into mode functions Uki(r, <p), which are 
eigensolutions of V^zifcj(r, ip) = 0 , 


£(r ,t) = Yu kt i(r, <p)£k,i CM)- (5) 

k,l 


The mode functions 


Uki(r, ip) 


Cu 

Wo 



/■WQ+iktp j^\k\ 


Wo 



, ( 6 ) 


are complete orthogonal set in (r, <fi). Lf are the associated 
Laguerre polynomials and Cki are appropriate normalization 
constants. I and k denote the radial and azimuthal index of 
the mode functions, respectively. The decomposition ([5]) is 
adequate as it describes Gauss-Laguerre modes of paraxial 
light propagation l23l (where wq wq (z)) for z values well 
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within the Rayleigh length zr = j- w'q of the focal plane, 
where wq(z) ~ wq = const. 

We can decompose the optical polarizations a ge (r,t) and 
<7 gr (r,i) into Uk,i(r, ip) in an analogous way and obtain from 
(|4} and the completeness of the mode functions 


where we defined the effective decay rate of the bright-state 
polaritons as r e ff = fl 2 /r, with fig = g 2 n + fi 2 . Comparing 
the off-diagonal coupling between dark- and bright-state po¬ 
laritons to the difference of the diagonal terms, we find that we 
can treat the off-diagonal terms perturbatively if the condition 


(d t + cd z ^j£ki(z,t) = 

[ drd vV n ( r ) u ki( r > <P)u mn (r, ip)a™’ n (z, t). (7) 

m,n 

If the atomic density n( r) is slowly varying spatially in r on 
the scale wq and furthermore is independent on ip, orthogo¬ 
nality of the modes yields 


c\k\ < |r e fj| 


(ii) 


is fulfilled for all relevant fc-values of the polariton field. We 
note that this condition sets a limit for the characteristic length 
scale of the dark-state polaritons 


W » § « ^iabs, 
7 


( 12 ) 


(d t + cdPj £ kl (, z,t) = igy/na^ ( z,t). (8) 

If furthermore the driving-field Rabi frequency is independent 
on and slowly varying in r, the Heisenberg-Langevin equa¬ 
tions ([3]> decouple in the transverse modes u ky i{r , <p). In this 
case eq. 0 can be reduced to a one-dimensional propagation 
equation for individual transverse modes. Note that this is 
only correct as long as interactions are disregarded. The ef¬ 
fect of the latter will be discussed later on. We are interested 
in particular in the input mode £ 0 0 for which we drop the in¬ 
dex ( 0 , 0 ) for simplicity in the following. 


C. Slow-light polaritons 


Assuming that the control field f2(z, t) is slowly varying in 
z on the scale w o yields a system of linear partial differential 
equations for the probe field and the atomic coherences, the 
Maxwell-Bloch equations, which separate in the transverse 
modes. Restricting to the lowest transverse mode and defining 
x = (£, CT gr , <r ge ) T we can write these equations as 

( —icd z 0 —gy/n\ 

0 _ 0 -fi . (9) 

—g^/n —fi —*T J 

Transforming to Fourier space we can find the eigenmodes of 
the Maxwell-Bloch equations. Restricting ourselves to small 
momentum (k « 0 ) the eigenmodes are the so called dark- 
state polariton mode and two corresponding bright-state po¬ 
lariton modes 12T1 . The dark-state polariton is a superposition 
of electric field and atomic coherence tf gr according to the 
basis rotation 'k = cos 9£ — sin 9s/na S[ , where the mixing 

angle 9 is given by tan 2 9 = g 2 n/fl 2 . We define the bright- 

state polariton as <k = sin 6 >£ + cos 9y/na gr , which is not an 
eigenmode of ([9]), but is a more convenient definition. As¬ 
suming the complex decay rate F of the intermediate level to 
define the shortest timescale as |r| _1 , the optical polarization 
tfge can be adiabatically eliminated. Rotating the remaining 
fields £, <j gr to the basis y = (d7, <h) T leads to the equation of 
motion 


id t y = H'y, H' 


f c cos 2 9d z c sin 9 cos 9d z \ 

\^c sin 9 cos 9d z c sin 2 9d z + r e ff J ’ 

( 10 ) 


where the second approximation holds in the off-resonant 
case, and L a b s = g 2 n/c r ) denotes the resonant absorption 
length of the medium. It is immediately clear that many-body 
effects can only be observed in media with large optical depth 
OD= L/L a bs, with L being the medium length. 

From the equation of motion we can read off the effective 
Hamiltonian in second quantization 


H n = — ic / d z 


cos 2 9^>\z)d z ^!(z) + sin 2 9&(z)d z <h(z) 


+ 


sin 9 cos e(p\z)d z $(z) + $ t (z)a 2 4-(z)) 


+ —&(z)$(z) 
c 


(13) 


describing the free time evolution of dark- and bright-state po¬ 
laritons, respectively, and their mutual coupling. 


D. Rydberg interactions 

For highly excited Rydberg states | r), dipole-dipole inter¬ 
actions between atoms become important, due to their large 
dipole moments CD. Atoms in a state |r) interact with the 
van der Waals interaction potential V(r) = C§/\r\ 6 with in¬ 
teraction strength Cq. For an ensemble of Rydberg atoms the 
microscopic Hamiltonian describing the interaction reads 

V = \ & ™ V ( r i - T,)^, (14) 

».jV* 

A (i) 

where (T rr denotes the projection operator to the Rydberg 
state of atom i at position r,. The interaction potential © 
leads to a space-dependent level shift which is the dominating 
energy scale on small length scales. If the atoms sit too close 
to each other the transition from a single Rydberg excitation 
to the doubly-excited state by laser light is prohibited, leading 
to a strong suppression of a second excitation, which is the 
so-called Rydberg blockade ll25l . 

Assuming a small excitation probability per atom to the Ry¬ 
dberg state, which allows to set o rl « < 7 f j, r <fg r , transforming to 
coarse-grained operators according to equation <[T]) and per¬ 
forming the continuum limit as above leads to the continuous 
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interaction Hamiltonian 


r, r' and ip, ip', 


Hint = y Jd 3 rJd 3 r' V r (r-r , )o-| r (r)CTt r (r , )CT gr (r , )CT g r(r). 

(15) 

The atom-light interaction, eqs. <0 drives the atoms into 
stationary dark states, i.e., states from which there are no 
spontaneous emission losses, provided that the two-photon 
detuning 5 is sufficiently small. As a consequence the atomic 
medium becomes transparent (EIT), and light propagates un¬ 
depleted with reduced group velocity. The van der Waals 
interaction between Rydberg excitations gives rise to a two- 
photon level shift, which exceeds the EIT linewidth when the 
distance becomes less than the EIT blockade radius 


:=C, 6 Jdp J dp' j r dr j r'dr' 

0 0 0 0 

tt fc 1 l 1 ( r ) u fc 2 J 2 ( r, )“fc3J3(r , )«fc4J 4 (r) 

V 1 3 ' I 1 ”-* 

r 2 + r ' 2 + 2 rr' cos (ip — ip') + (z — z ') 2 


The angular integrals can be calculated analytically by residue 
integration yielding 


-trk 1 k 2 k 3 fe 4 

V hhhU 


(z 


z') 


■>k 1 — k i ,k 3 — k 2 


( 20 ) 


a B = y/|C 6 r|/H 2 . (16) 

As a consequence for short distances the mixing between 
dark- and bright-state polaritons becomes strong and the po- 
lariton picture is no longer adequate ESI E27). However, if the 
excitation density is sufficiently smaller than a ~^ 3 the polari- 
ton picture is expected to hold true. It should be noted that 
under slow-light conditions O 2 becomes small and thus the 
EIT blockade radius becomes large. In fact in the limit of 
light storage, where f 2 (i) — > 0 one finds a B — > oo and thus it 
is unclear if light storage in a Rydberg EIT medium is possible 
at all. We will show below, however, that the critical distance 
a c between excitations at which a significant mixing between 
polaritons sets in is not given by ( p~ 6 | ) but by 

a c = ^|c 6 r|/n 2 , (17) 

which stays finite in the limit of light storage, as fi 2 —> g 2 n. 


E. Reduction to a one-dimensional model 


Rydberg interactions can lead to a scattering between dif¬ 
ferent transverse Laguerre-Gaussian modes u^jir. p) and 
thus even the paraxial propagation of slow-light polaritons in 
a Rydberg gas becomes in general a three-dimensional prob¬ 
lem. In the following we will show that also in the presence of 
interactions the effective one-dimensional description is valid, 
provided that the beam waist of the beam is small compared 
to the Rydberg blockade radius a B . 

To derive the interaction Hamiltonian of an effective one¬ 
dimensional model we use the decomposition of the spin- 
flip operators <x gr (r) into Lague rre-G aussian modes <j gr (r) = 
; U/-i(r)a^(z) as in section IIB Thus we can decompose 
the interaction into different potentials describing interaction 
between different transversal modes 


H 


int 


= y/ dz / dz ' £ 


fcl &2 


x a^ lh (z)a^ 2h (z')a^ h (zf)a£ u {z), (18) 

where the effective potentials are defined by integrating out 


The further evaluation has to be done numerically. We 
are interested in the scattering processes of an initial Gaus¬ 
sian mode with zero angular momentum into higher order 
Laguerre-Gaussian modes which are governed by the poten¬ 
tials V^“ fc 0 00 . The modes uu with k ^ 0, i.e., higher-order 
azimuthal modes, have vanishing amplitude at r = 0 and the 
corresponding interaction processes are suppressed compared 
to the k = 0 modes. Thus we will restrict the following 
discussion to scattering processes into azimuthally symmet¬ 
ric modes. In Fig. [2] we show different potential curves for 



Figure 2. Different interaction potentials between two pho¬ 

tons initially in Laguerre-Gaussian modes Z 3 , Z 4 sitting at positions 
( 2 , z' = 0) and finally in Laguerre-Gaussian modes l\, I 2 , as indi¬ 
cated by the inset. In particular we choose initial Gaussian modes 
l 3 = Z 4 = 0. The dashed lines indicate power-law fits in the 
regime wo < z < zr obtained by assuming zero beam divergence 
w(z) = Wo- 


the interaction of two photons at relative distance z, initially 
in the Gaussian mode I 3 = I 4 = 0 and finally in (higher or¬ 
der) Laguerre-Gaussian modes with radial indices ( 1 ,( 2 - The 
relevant reference is the forward scattering potential VJ^o'oo 0 - 
We find that for distances smaller than the beam waist wq 
all potentials show a z ~ 4 dependence with different ampli¬ 
tudes due to the small overlap of the different modes. At 
z = wo the potentials cross over to a power law z~ a , where 
we find numerically that a ~ 6 + 2 (Zi + Z 2 ), before they 
flatten again for z > z B . For strong interactions small dis¬ 
tances are blockaded, where typically the blockade distance 
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is on the order of a few wq. Therefore the relevant regime 
is given by wq < as < z < zr, where the modes func¬ 
tions Uki{ r) are approximately constant in ^-direction and the 
potentials scale as z~ a . One recognizes that already for dis¬ 
tances slightly larger than the beam waist the scattering am¬ 
plitudes into higher Laguerre-Gaussian modes are orders of 
magnitude smaller than the forward scattering one. (Note that 
the potentials describing scattering between different excited 
modes Uki are suppressed by at least an order of magnitude 
more and decay much faster compared to the interaction be¬ 
tween the initial Gaussian modes.) Thus we can safely neglect 
the scattering processes into higher modes if the distance be¬ 
tween excitations is sufficiently larger than the blockade ra¬ 
dius as > wq. This yields an effective one-dimensional in¬ 
teraction described by the van der Waals potential. Omitting 
again the indices 0 we find 

Hint = y JdzJdz'V{z- z')ol r (z)al r (z')a gr (z')a gr {z) 

( 21 ) 

Note that the ID operators <r gr (z) = <7™(z) have a different 
physical dimension as the spin flip operators cj gr (r) in three 
dimensions. 


F. Effective model for Rydberg dark-state polaritons 


reservoir in the vacuum steady state. The free time evolution 
of the polaritons in zeroth order in the coupling is described 
by the differential equations 

d t ^(z,t) = — ccos 2 8d z d>(z,t) 

d t $(z,t) = —csin 2 6d z <f>(z,1) - T eff $(z,f), 


where we have disregarded the Langevin noise terms. As can 
easily be derived, the solution of these equations is given by 
^(z, t) = ^(z — ccos 2 8t, t— t) and $(z, t) = e~ reffT< l>(z — 
csin 2 6t, t t). Thus the free dark-state polaritons propagate 
with group velocity v g = c cos 2 8 with a stable shape, while 
the bright-state polaritons propagate with c sin 2 8 and are sub¬ 
ject to decay with the rate Re[T e ff]. Assuming that |T e ff | —1 
defines the fastest time scale, justifies adiabatic elimination of 
the bright-state polaritons. If there is no external driving of 
bright-state polaritons, the steady state is the vacuum and we 
find that bright-state polaritons are delta-correlated in space 
and their correlation functions decay in time as 


($(x,t)&(y,t- t)) « e r ° ltT 6(x-y), r > 0, (23) 
(I ,t)& (y,t - r)4> t (y,f - r)) « 


2r 0 ffr 


S(x — y — vt)S(x' — y' — vt) 


+S(x' — y — vt)S(x — y' — vt) 


, t > 0. (24) 


Transforming the ID interaction Hamiltonian ( |2T| to the 
polariton basis according to y/nd gI {z) = —sin8 d> (z) + 
cos $<l>(z) yields 

H int = i/ dzj dz , y(z-z , ){[sin6»^ t (z)-cos6»$ t (z)] 

x [sin d^(z') - cos6>4> t (z')] [sinfl^z') - cos6><l(z')] 

x [sin (z) - cos 6>$(z)] j. (22) 

After expanding one can identify terms describing: (i) a two- 
body interaction of dark-state polaritons with relative strength 
sin 4 8, (ii) an interaction of bright-state polaritons with rela¬ 
tive strength cos 4 8 and (iii) nonlinear coupling terms of dark- 
and bright-state polaritons. Note that in the slow light regime 
1 « sin 4 8 2> cos 4 8 the dark-state polaritons, consisting 
mainly of Rydberg excitation, are strongly interacting while 
the bright-state polaritons, consisting mainly of electric field 
excitation, are weakly interacting. In lowest order in cos 0 
only a nonlinear interaction term between Rydberg dark-state 
polaritons remains. In the following we want to consider also 
the dominant correction terms, caused by the coupling of the 
dark-state polaritons to bright-state polaritons, which we can 
assume to be in a vacuum state. 

The full time evolution of dark- and bright-state polaritons 
is governed by the sum of free and interaction Hamiltonian 
and can analogously be decomposed into dark-state, bright- 
state and coupling terms, H = Ho + Hint = H^ + H$ + H^$. 
We make a standard system-plus-reservoir approach |[28ll and 
derive an effective Liouvillian for the Rydberg dark-state po¬ 
laritons by treating the bright-state polaritons as a Markovian 


All remaining correlation functions of bright-polariton opera¬ 
tors vanish in the vacuum state. 

In an interaction picture with respect to the system and 
reservoir Hamiltonian I I ,j, + H$ the dark-state polariton de¬ 
grees of freedom are described by the density matrix p = 
tr$x, which we get by partial trace of the full density ma¬ 
trix x over the bright-state polaritons. The time evolution of p 
is then governed in Born approximation by the equation 


Pv(t) = - J dr tr$ | [H*$(i),p(t) <g> /£>$]] j , 


(25) 

where p,\, denotes the bright-state polariton steady state den¬ 
sity matrix. We note that the free part of 11corresponds to 
a transformation to a frame comoving with the group velocity 
v g = ccos 2 8. Introducing the operator 


L = — sin 4 8 cos 8 


[ dsV r (z—s)^^(s)'T(s )+*—y 
J sin 8 


'T(z), 


(26) 

allows us to write the system-reservoir-coupling Hamiltonian 
in the interaction picture in the compact form 


H^$(i) = jdz{&(z)L(z) + m(z) + Jdz'V(z — z') 

x sin 2 8 cos 2 8 ^4 t (z,f)4> t (z , ,<)4'(z , ,f)T'(z,f) 

+ }, (27) 

where we omitted terms third order in the bright-state polari¬ 
tons since the three particle correlations of bright-state polari¬ 
tons are vanishing in the vacuum. We insert Hamiltonian (|27j) 
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into equation ( |25j ) and follow the standard derivation of a Mas¬ 
ter equation 1281 . By performing the Markov approximation 
we get a Master equation in Lindblad form describing an ef¬ 


fective time evolution of the dark-state polaritons. After trans¬ 
forming back to the moving frame Schrodinger picture we ar¬ 
rive at 


p=i K r z 


p,i)(z)L(z) +i dz dz'V(z-z') p, &(z)'fr 1 '(z')'I'(z')\l/(z) 


. A sin 4 9 cos 4 9 


fl? 


+ 


dz Jdz'V 2 (z-z r ) p, d' t (2)'J' t (2 / )#(2; , )'k(2) + ^ jdz (2L(z)pL\z) - jp, D{z)L(z)^j 

r dz f cl z' V 2 (z - z') (2&(z')V(z)p&{z)&(z') - {p, (*')$(*')*(*)}) (28) 


7 sin 4 9 cos 4 9 


where {-,-} denotes the anti-commutator. The coupling of 
dark- and bright-state polaritons gives rise to various effective 
unitary and dissipative terms in equation ( |28j ) influencing the 
evolution of dark-state polaritons. Specifically, we find the 
following unitary terms: a kinetic energy term with an effec¬ 
tive mass defined by 

m -1 = 2u g cA ^ = 2u g —T a b S sin 4 9, (29) 

7 

a drift term ~ f dz'V(z — z')& (z')^(z') - mediated by in¬ 
teractions - giving corrections to the group velocity v g , and, 
finally, higher order terms in the interaction, namely a three 
body interaction and corrections to the two-body interaction 

~ V 2 (z — z'). 

Moreover, since the bright-state polaritons decay the cou¬ 
pling also leads to effective loss channels for the dark-state 
polaritons in equation ( |28| . The decay processes are described 
by the operator L defined in ( |26| . This operator consists of 
two terms arising from the two effective loss channels de¬ 
scribing coupling between dark- and bright-state polaritons, 
namely a linear coupling outside the EIT window and a non¬ 
linear coupling arising from the interaction. Both terms lead 
to losses of dark-state polaritons. In particular, we identify 
a generalized single particle loss generated by the Lindblad 

operator r/r(z)^(z) with the operator valued loss rate 


f(z) = 


7 cos 2 0 sin 6 9 


dx V(z — tc)^ , ^(a;)4 r (a;) 


(30) 


Note that this loss as well as the corrections to the drift term 
cannot be simply treated in a mean-field approximation, since 
V(z — z') can assume arbitrary values. However, the full 
Lindblad operator gives rise to a space-dependent loss process 
~ r -12 that is strong for small distances r. 


G. Effective Hamiltonian 

After an initial transient we expect polariton states to evolve 
according to equation (|28| and to show an blockaded volume 


for small distances. We assume that the product of interac¬ 
tion energy and mass is always positive, corresponding to re¬ 
pulsive interaction that prevents excitations to propagate into 
the blockade volume. On large length scales the interaction 
decays as V(r) ~ |r| -6 , thus the effective interaction terms 
linear in V are the dominant contributions and we may ne¬ 
glect the two body interaction with potential ~ V 2 ~ M -12 
and the three-body interaction as well as the loss term with 
rate ( |30| i. As a second approximation we assume that we are 
in the slow-light regime, where cos <C 1 and only excita¬ 
tion energies inside the EIT window are allowed. Finally, we 
consider an off-resonant driving scheme where the fields are 
far detuned and 7 <C A which allows us to neglect the re¬ 
maining loss terms. An expansion in cos 9 up to second order 
reduces the master equation to a von Neumann equation for 
the density matrix, governed by the Hamiltonian 


H = - Jdz&{z)^-$(z) 


+ Cq sin 4 9 dz dz 




Z - 2 ' 


v|6 


(31) 


in the moving frame Schrodinger picture. 


III. TWO-PARTICLE SIMULATIONS AND 
VERIFICATION OF EFFECTIVE MODEL 


We are interested in the effect of the Rydberg interactions 
on the evolution of propagating dark-state polaritons. For 
two copropagating polaritons this has been considered theo¬ 
retically and simulated numerically in l29l and, in a resonant 
as well as off-resonant regime recently been demonstrated in 
experiments mum. We want to extend this to the strongly 
interacting far off-resonant regime, where we predicted in the 
previous section that the dynamics is governed by the effective 
Hamiltonian ( [3T| in a comoving frame. 

It was shown in ||29l that the van der Waals interaction 
leads to an avoided volume for the propagating photons of size 
2a b- Near single-photon resonance, | A| < 7 , two photons at 
distance smaller than get absorbed and an incident two- 
photon wave-packet evolves into a non-classical state that is a 
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Figure 3. Splitting of initial Gaussian wave packet due to van der Waals interaction. Full two-excitation wave function simulation for the 
parameters g = lOfl = IO 7 and A = 47 . The optical depth per blockade is chosen ODb = ob/L abs = 7. The color plots in a) and b) 
show the two-excitation probe field amplitude EE(z\, z 2 ,t) at times t = 0 and t = 0.5 L/v s corresponding to the times where the center of 
the incident Gaussian wave packet is at the medium boundary and at the middle of the medium, resp. Plot c) shows a cross section (blue solid 
line) of the middle color plot along the diagonal r = z\ — z 2 for 2 R = zi + z 2 =L = 100 Ob, as indicated by the dashed-dotted line in b). We 
compare the cross section to a pulse propagating in absence of interactions (Gaussian, purple dotted line) and the two-particle ground state of 
the effective Flamiltonian ( |3 1 [ > (yellow dashed line),where we impose periodic boundary conditions with a period corresponding to the distance 
of the two peaks (blue solid line). The red dash-dotted line shows a thermal two-particle state of the effective Flamiltonian. For illustration 
purposes we continued the two-particle ground state and the thermal state. 


statistical mixture of a single excitation and a correlated train 
of (two) photons separated by 2a b 1291 . In the off-resonant 
regime (|A| 7 ) the absorption plays no role but two pho¬ 

tons closer than the blockade distance propagate with the vac¬ 
uum speed of light and thus escape. Moreover, the repulsive 
interaction prevents the photons to get inside the blockade ra¬ 
dius. 

To validate the simple picture we perform numerical sim¬ 
ulations of the full ID Maxwell-Bloch equations for two co¬ 
propagating polaritons subject to van der Waals interaction. 
The two-excitation wave function | SI /2 (^)) is composed of 
six two-excitation modes describing all combinations of the 
single-excitation modes E, a gI and er ge . Note that in presence 
of decay the norm of the wave-function is not conserved. 

We consider an initial Gaussian two-photon wave packet 
(l^o) = £2£T |0)) in free space, where we assume that the 
spectral width fits well inside the EIT window of the medium 
in absence of interactions. The state vector at a later time t is 
then described by two-excitation wave functions 

|^(*)) = fdzi jdz 2 Y, FG(z 1 ,z 2 , (zijG^Q^lO) 

F,G 

where FG £ {EE, PP, SS, EP, PE, ES, SE, PS, SP} de¬ 
note all possible two-particle wave functions corresponding 
to excitations in the electric field (E), the optical polariza¬ 
tion (P), and the spin polarization ( S ). We then simulate the 
propagation of the pulse from free space into the Rydberg gas 
by considering a space-dependent atomic coupling strength 
g 2 n = g 2 n(z), where we use a step function for n(z). To take 
care of the EIT pulse compression for small group velocities 
we use an adaptive spatial grid spacing. In Fig. [3]a), b) we 
show two snapshots of the time evolution of the EE(z\, z 2 )~ 
component incident into the medium at the lower left corner 
(zi /2 = 0) and propagating along the diagonal to the upper 


right corner of the pictures. As can be seen from the figure 
in the considered off-resonant regime the wave packet splits 
at the medium boundary into two parts separated by the off- 
resonant blockade radius as, indicating a photon antibunch¬ 
ing as shown for the resonant case in l29l . However, af¬ 
ter longer propagation inside the medium the separated parts 
move even further apart due to the repulsive interaction on 
length scales Az as- 

The results in Fig.[3]verify the expected time evolution gov¬ 
erned by Hamiltonian ( |3T| ). This can be seen from Fig. [3]c), 
where we compare a cross section (blue solid line) of the 
evolved pulse to the two-excitation ground state of the effec¬ 
tive Hamiltonian (yellow dashed line) calculated by employ¬ 
ing periodic boundary conditions. For illustrative purposes we 
have periodically continued the ground-state wave-function in 
the plot (yellow dash-dotted line). 

We recognize two important things: First, we find a re¬ 
markably good agreement of the numerically evolved wave- 
packet to the ground state of the effective Hamiltonian up to 
a length scale of ~ ±25 as- Secondly, contrary to Ifl2l . we 
do not find a localized two-particle excitation around r = 0 
(bound state). That the final state of the evolution is so close 
to the ground state is somewhat unexpected since the initial 
wave-packet contains many highly excited components. The 
presence of a finite EIT frequency window sets some upper 
bounds on the energy of the polaritons. Since the system is 
not integrable, we can assume that the polaritons quickly ther- 
malize after entering the interacting medium. However, as can 
be seen from comparison with the red (dash-) dotted curves 
in Fig. [3] estimating the temperature by the EIT transmission 
window HD is much too high. A possible explanation of this 
is the presence of a mechanism similar to evaporative cooling. 
The high-energy modes of the Rydberg dark-polariton modes 
are preferentially coupled to bright-polariton modes and de¬ 
cay or propagate away. A full analysis of this mechanism is 
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however not subject of the present paper. 

Let us briefly comment on the fact that we do not see 
bunching of photons in the numerical simulations of the full 
Maxwell-Bloch equations as in HU- We here consider the 
case mCe > 0, where the effective Hamiltonian ( |3 11 does not 
have bound states which could explain the observed bunching. 
The full set of Maxwell-Bloch equations on the other hand 
should reproduce the bunching. However, our simulation is 
run in a regime with much larger optical depth per blockade 
ODb as in m. As will be discussed elsewhere l27l . in this 
regime not a single, but many bound states exist which cannot 
be excited, however, by the incident photon pulse. For values 
of ODb much smaller than considered in the present paper we 
also find bunching in our numerics. 

IV. GROUND STATE PROPERTIES OF EFFECTIVE 
HAMILTONIAN AND LUTTINGER LIQUID APPROACH 

We have seen in the previous section that a finite-length 
light pulse entering an EIT medium with Rydberg interactions 
with large optical depth per blockade ODb generates a many- 
body state of Rydberg dark-state polaritons with an energy 
close to the ground state of the effective Hamiltonian ( |3T| in 
the moving frame. In this section we will discuss the ground- 
state properties of this Hamiltonian using numerical DMRG 
simulations as well as a Luttinger-liquid approach. 

A. DMRG ground state simulations 

The physics of © is governed by an interplay of kinetic 
energy contribution trying to delocalize the massive particles 
and polariton-polariton interaction leading to spatial order. 
Note that the effective mass ( |29| ) can be negative, depend¬ 
ing on the sign of the detuning A. We assume in the fol¬ 
lowing that the product of mCe has always a positive sign, 
corresponding to a repulsive interaction. To analyze the inter¬ 
play between kinetic and interaction energies we calculated 
ground state correlations using the density matrix renormal¬ 
ization group l30l . Specifically we calculate single- and two- 
body correlation functions for different interaction strengths, 
quantified by the dimensionless parameter 0, which is defined 
as the ratio of interaction energy at average inter-particle dis¬ 
tance 1 1pQ compared to the Fermi energy fl40l . 

0= (A^(7) 2 OD ^ (32) 

In Fig. [4] we show density-density and first order correlations 
calculated by DMRG. For short distances the density-density 
correlations are strongly suppressed, i.e., g^ 2 \0 ) = 0, indi¬ 
cating the photon blockade of two copropagating excitations. 
For small 0 corresponding to weak interactions the first order 
correlation functions govern the long-range behavior, indicat¬ 
ing a superfluid state. Increasing the interaction strength leads 
to a strongly pronounced and slowly decaying density wave 
(CDW) with fast decaying first order correlations. From the 
plot we can extract a power law decay z - ;i for distances larger 



Figure 4. DMRG correlation functions. Top: 1-density-density cor¬ 
relation functions for different interaction strength in a double loga¬ 
rithmic plot and density-density correlation functions in a linear plot 
in the inset. Taken from CD- Bottom: corresponding first order cor¬ 
relation functions. 


1/po- Note that for zpo > 4 finite size effects influence the 
results. Due to the low dimensionality (d = 1) of the model 
no slower-than-power law decay of correlations can be found, 
i.e., the model cannot exhibit tme crystalline order. However, 
the latter can in principle be created by engineering an addi¬ 
tional lattice potential, i.e., a space-dependent two-photon de¬ 
tuning for the polaritons leading to a sine-Gordon like model 
for commensurate fillings. This model exhibits a quantum 
phase transition to a gapped phase with true crystalline or¬ 
der GU [32 - 


B. Luttinger liquid approach 

The low energy physics of a gapless one-dimensional sys¬ 
tem can be described in terms of a Luttinger liquid (LL) 
model. Assuming that we have a state of N polariton exci¬ 
tations in a frame moving with group velocity v s distributed 
over length L yields an excitation density of p 0 = N/L. 
Working in the far off-resonant regime, spontaneous decay 
can be neglected after an initial transient and the repulsive in¬ 
teraction suppresses the coupling to the bright-state polariton 
loss channels. Therefore the number of polaritons can be as¬ 
sumed to be constant. Following standard LL theory ll32l we 
express the bosonic polariton field 'l'(z) by fields <!>{z) and 
6 {z) representing density and phase fluctuations of the long- 
lived sound-modes, respectively. The density of the polariton 
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field transforms as p(z) = 'K (z)’T(z) = —d z (j)(z)/ tt. Writ¬ 
ing the system Hamiltonian in terms of and the canoni¬ 
cally conjugate mode n(.z) = d z 0(z)/ir yields the LL Hamil¬ 
tonian 


HLL= 2ir z 


v,K(*n(z)f + 


(33) 


The Hamiltonian is fully defined by two parameters, the speed 
of sound v s and the Luttinger parameter I\. The K param¬ 
eter universally determines the long-range behavior of the 
ground state correlation functions. In particular, the equal 
time density-density correlations are given by 132) 


(pO)p(0)) = Po - ^2 ~2 + A iPo cos(2np 0 z)e 2G , 

(34) 


where G^iz) = ( \(j>(z ) — </(0)] 2 ) and A2 is a non-universal 
amplitude that has to be determined numerically from the mi¬ 
croscopic model. We are particularly interested in the last 
term that oscillates spatially with the period 1 /po correspond¬ 
ing to a charge-density wave (CDW). The spatial decay of 
the CDW is governed by the correlation function (-r 2C,i, A z ) 
which in the ground state is given by a power law (a/z) 2K 
with exponent 21 \ and a short distance cutoff a, which we 
choose to be the shortest length scale which is typically the 
average distance of excitations p^ 1 . In contrast to this the 
ground state first order correlation function describing super¬ 
fluid order is given by 


($t( 2r )$( 0 ))=p 0 A 1 (^J , (35) 

where the amplitude A\ is also non-universal. 


C. Quasi-Crystalline state under stationary EIT conditions 

The long-range behavior of first order ( |35| i and density- 
density correlations are in the ground state both given 
by a power law with exponents 1/2 AT and 2 K, respectively. 
Comparing this with Fig. [4] shows that they both agree qual¬ 
itatively well with the DMRG results. Depending on K ei¬ 
ther superfluid (first-order) (K 1/2) correlations or the 

CDW (K -C 1/2) dominate the long-range nature of corre¬ 
lations, where the point K = 1/2 marks the crossover point, 
where both correlations decay with an exponent 1. For real¬ 
ization of a quasi-crystalline state a regime A' C 1 has to 
be realized. Notice that the long-range interaction yields a 
A'-parameter smaller 1 for a sufficient interaction strength, 
indicating more pronounced density correlations than free 
fermions which is the strongest that can be achieved for 6 - 
interacting bosons l32l . e.g., photons interacting via a Kerr 
nonlinearity ll33l 1341. 

For our model ( |3T| the dependence of the A'-parameter as 
a function of microscopic parameters cannot be given analyti¬ 
cally, but has to be determined numerically. We can extract A' 
from DMRG simulations, where we calculated the compress¬ 
ibility x _1 = pl§f~ Q = P 2 0 L§£. From this we get the ratio 
K/v s = 7 tpqX and together with the product v s K = npo/m 



Figure 5. A'-parameter as a function of the optical depth per blockade 
ODb for A /7 = 10. Solid curves are interpolated numerical data 
from DMRG results, corresponding dashed curves are approximate 
values according to m The three different colored pairs of curves 
are for different values podB as indicated in the plot. The dashed 
vertical line indicates experimental values from CD. 


which is constant for any Galilean invariant system lf32l we 
can extract the A'-parameter as a function of the dimension¬ 
less parameter 0, ( f32] >. In the case of power law interactions 
an approximate analytical formula for the AT-parameter has 
been given in l35l . which is asymptotically correct for small 
and large 0, 


K = 




- 1/2 


(36) 


From this we can estimate the experimental requirements 
for reaching the strongly interacting regime (I\ <C 1) under 
stationary EIT conditions. An important restriction of station¬ 
ary EIT is set by the Rydberg blockade: When two excitations 
get closer than the EIT blockade radius as, eq. 0. they are 
either absorbed or transformed into fast propagating bright- 
state polaritons which escape. Thus the excitation density is 
limited due to photon blockade to values 


Po^b < 1 


(37) 


and in the regime of small K we have K ~ OD/ 1 indicating 
that even for a pqcib ~ 1 the required optical depth is orders of 
magnitude higher than experimentally feasible imcf.Fig.0 
Moreover, increasing OI) I; by changing the blockade distance 
as requires a smaller excitation density p () for fixed an d 
thus limits the number of excitations in a medium with a finite 
length such that in the limit ODb A> 1 only a single excita¬ 
tions would be allowed inside the medium for stationary EIT 
driving. 








V. REACHING THE STRONGLY INTERACTING 
REGIME: LIGHT STORAGE IN RYDBERG GASES 
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A. frequency pulling in adiabatic slow-down of Rydberg 
polaritons 

The problems to reach the required conditions for quasi¬ 
crystalline order of Rydberg polaritons under stationary EIT 
conditions can be overcome using a dynamical protocol, i.e., 
by considering light storage or dynamical slow down of po¬ 
laritons while inside the medium. One recognizes from © 
that during slow down and ultimately light storage the EIT 
blockade radius ~ fi -1 / :i would diverge. Naively one 
would expect that as a consequence the smallest possible dis¬ 
tance between Rydberg polaritons diverges as well, when the 
group velocity goes to zero. Remarkably this is not the case. 
To see this let us consider the dynamics of a dark-state po- 
lariton during storage with a finite two-photon detuning Sq, 
resulting, e.g., from a second nearby Rydberg excitation. As 
has been shown in j36l a small two-photon detuning causes a 
time-dependent phase shift (chirp) of the dark-state polariton 
during slow down: 


4 /(z,t) = T — c J drcos 2 (r), 0 ^ 

x exp |zi5o J drsin 2 (r)j>. (38) 


As a consequence the spectrum of the probe field £(z,t) = 
cos 0(t)'Sf(z,t) assuming a slowly varying mixing angle 9(t) 
can be expressed as 


S(z,u)) = J dre lWT (£\z,t)£(z,t — r)^ 

— OO 

= (39) 

One recognizes two effects: First, there is a spectral narrow¬ 
ing proportional to cos 2 9{t), which guarantees that during 
light storage the pulse spectral width remains less than the 
EIT transparency window if it did so at the beginning of 
the storage process EH. Secondly, and more importantly in 
the present context, there is a pulling of the center frequency 
of the pulse towards the two-photon resonance 

S(t) = 5o cos 2 9{t). (40) 

This effect has been observed experimentally in fl38ll and is 
responsible for the fact that the two-photon linewidth of EIT 
light storage ll36l 



Figure 6. Illustration of transmission spectrum T(w) of the EIT 
medium (top) and spectrum of the electromagnetic field ( bottom ) in 
arbitrary units as function of the EIT mixing angle 6 during a dynam¬ 
ical light storage protocol. 


As a consequence of this the minimum distance of two 
Rydberg polaritons is not given by the EIT blockade dis¬ 
tance, ( p~6l >, but by the critical distance 

a c = V\c 6 r\/ni (42) 


We have verified this by two-particle simulations, which are 
illustrated in Fig. [7] Shown is the spin-spin component of an 
initial two-photon wave-packet propagating in an EIT medium 
with Rydberg interactions as in section III After an initial 
transient an avoided volume is formed corresponding to the 
blockade radius o^(0 0 ) given by the initial drive-field Rabi- 
frequency flo and larger than the critical distance a c . When 
the drive field is adiabatically turned to zero as shown in the 
inset, the EIT blockade radius as (12(f)) increases and even¬ 
tually diverges. The avoided volume stays however approxi¬ 
mately constant. 

Therefore adiabatic slow down of Rydberg polaritons al¬ 
lows to increase their effective mass and thus the ratio of in¬ 
teraction energy to kinetic energy, 0 , without reducing the 
polariton density pq. In this way it is possible to enter the in¬ 
teresting regime of strong interactions between Rydberg dark- 
state polaritons! 


B. Time-dependent Luttinger liquid approach 


S 2ph = n 2 J |T| (41) 

is determined by the collective Rabi frequency fl e rather than 
the control-field Rabi frequency (Here T = 7 + i A and 
we have set the optical depth OD = 1.) 


In the previous section we have shown that by employing a 
time-dependent protocol a regime of strongly interacting po¬ 
laritons with I\ ■Cl can be reached, as opposed to a contin¬ 
uous driving. Through two-photon simulations we have val¬ 
idated that during storage the polaritons are separated by the 
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Figure 7. Cross section of S'5-component during storage. The con¬ 
trol field is turned off according to a protocol tanh(f/r), which is 
shown in the inset. The main figure shows the cross section of the 
doubly excited spin component for 3 different times, before starting 
the storage protocol, at time t = 2r , where 17(f) = |flo and at the 
end of the protocol, where Q(t) is approximately zero. 


finite distance a c . The ground state in the regime K <C 1 is 
strongly correlated as we have shown in section [TV A| How¬ 
ever, since the Luttinger liquid is gapless the quench from 
weakly to strongly interacting regime cannot be fully adia¬ 
batic. Nevertheless, as will be shown in the following a state 
with long-range CDW correlations that extend over a finite 
length can be prepared. 

To analyze the many-body polariton dynamics we calculate 
the time evolution of the correlation functions ( |34| starting 
from a given (weakly interacting) initial state which we as¬ 
sume to have been prepared under stationary EIT conditions. 
To get the correlation functions after storage we apply a time- 
dependent Luttinger theory as in ll39l . 

Storing the polaritons by tuning the control field SI to zero 
leads to time-dependent parameters K(t) and v s (t), where 
we can calculate their time-dependence from the effective 
mass analytically using ( [36] ) and also numerically from exact 
DMRG simulations (Fig.[5F. We follow standard bosonization 
procedure and decompose the canonically conjugate fields >1) 
and n into bosonic momentum modes [b q ,bj,\ = 5 PA . This 
transforms the Hamiltonian ( |33| i to 

w{t)b\b p - ^ ( b\b ] _ p + b_ p b p ) , 

(43) 

where w(t) and g(t) are given by the LL parameter as K{t) ± 
iT -1 (f), respectively. The time evolution of the system can 
then be calculated by solving the Heisenberg equations of mo¬ 
tion for the bosonic modes which are given as 


Hll = 


V s (t) 


\p\ 


id t 



VsttM (wit) —g(t)\ ( b p \ 

2 \g(t) -wit)) \bl p ) ' 


(44) 


This is a system of coupled differential equations each cou¬ 
pling the momentum modes b p and b] p . Performing a time- 
dependent Bogoliubov transformation b p = u P it)b P i 0) + 
Vpit)bL p i0 ) and an analog transformation for b'_ p maps the 


operator time dependencies to the Bogoliubov coefficients 
Up it) and v p (t). Applying these transformation to the Heisen¬ 
berg equations (|44| leads to coupled differential equations for 
the coefficients. Defining R p = (it p (f), v P it)Y these can be 
written as 


id t Rp = M P it)R p: Mp (t) = *®& ■ (45) 

In general these equations cannot be solved analytically. Since 
the matrix M p (t) is time-dependent, a transformation S p 
which diagonalizes M p (f) is in general time-dependent it¬ 
self and thus leads to nonadiabatic corrections. These correc¬ 
tions introduce again an off-diagonal coupling of the trans¬ 
formed variables = S ,_1 (f)R p . To find the correc¬ 

tions due to this off-diagonal coupling we calculate the ma¬ 
trix S p such that S p l M P it)S p = M^it) is diagonal and 
transform equation ( |45| accordingly. We find that the diag¬ 
onal entries of the resulting matrix M^it) are then given 
by ±|p|t> s (i) whereas the off-diagonal couplings S^Sp are 
given by —?'AT(f)/2A"(i). The latter can be neglected if 


\p\ » 


Kjt) 

2 Kit)v„(t) 


(46) 


For a fully adiabatic dynamics this condition has to be fulfilled 
for all relevant modes p at all times. However, this is impos¬ 
sible for all momentum modes since the Hamiltonian ( |43| is 
gapless. Nevertheless, bounding the timescale of changing 
Kit), i.e., of storing the polaritons, gives rise to a momentum 
scale p c such that all large momentum modes \p\ > p c ful¬ 
fill ( |46| ) at all times. We divide the set of momentum modes 
into two parts |p| > p c and \p\ < p c which gives rise to 
two spatial regimes. On small length scales corresponding 
to large momentum modes the correlations can adiabatically 
follow the change of parameters while on large length scales 
corresponding to small momentum modes they stay frozen in 
the initial state. This simple argument yields a length scale 
L 0 ~ 1/pc marking the crossover between the two regimes. 

A better estimate for the crossover-length between adia¬ 
batic and diabatic behavior can be found as follows: Corre¬ 
lations in real space propagate at the speed of sound v s , i.e., 
turning v s ~ yjI\ it) to zero during light storage freezes the 
correlations. Performing the storage in finite time lets corre¬ 
lations propagate up to a certain length scale L colr marking a 
crossover between an adiabatic regime, where the correlations 
follow the storage to a diabatic regime which shows the initial 
correlations. This length scale is given by the integral over the 
speed of sound 


d tv s (t), (47) 

which we will calculate for a particular time-dependence in 
the following. 
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C. Special solution 


D. Correlation functions 


Let us now discuss a particular time-dependence of the 
parameters, which allows for a semi-analytic solution of 
the time-dependent Luttinger model. To this end, we use 
the approximate dependence of K on the mass from equa¬ 
tion ( |36| to introduce a time-dependence of the AT-parameter 
via 0(t)/0 o = fl(0) 2 /fl(t) 2 . Furthermore, we iterate the di- 
agonalization procedure, i.e., we find a Matrix that diag¬ 
onalizes Mp . The off-diagonal corrections in this second- 
order adiabatic diagonalization are then given by the time 
derivative 


± 


d 
d t 



(48) 


For a Galilean invariant system we can express the speed of 
sound in terms of K and m via v a = irpo/mK. Plugging 
this into ( [48] ) and making use of the fact that AT and m have 
a unique relation, we can find a special time dependence of 
K[t ) (or equivalently m[t)) such that expression ( |48] > van¬ 
ishes for all times. This yields a differential equation for K(t) 
which has the solution 


A"(f) = e" acosh(t/r+c) , (49) 

where C = fAo + 1/ATo). With this we can calculate the 
critical momentum p c defined by the expression on the right 
of the inequality (|46|, 


Pc = 


Kit) 

K{t)v s {t) 


1 1 

-= — = const ., 

Trp 0 r]T L 0 


(50) 


which is now constant in time. We can express the crossover 
length L 0 = fg (poas) 5 v g (0)ODB |^yr in terms of the pa¬ 
rameters poa° B and OI) ( /> describing the initial conditions. 

For the maximal distance up to which correlations can build 
up during the storage protocol we find 

Acorr(i) = j Q ds V ^ S '> = Y 111 (j^j) ‘ (51) 

Apart from logarithmic corrections the maximal distance is 
given by Lq/2, allowing in principle for a large correlation 
length A C orr if AT(f) <C AT(0). However, for large times t the 
AT-parameter vanishes like A' (t) —> f i.e., it takes long times 
to reach small AT values. This limits the maximal A C orr(f) in 
practical realizations of the protocol. 

Switching off the control field such that AT gets tuned to 
zero according to ((49), the off-diagonal coupling of the twice 
transformed equations vanishes and we can integrate the dif¬ 
ferential equations to get the solution 


To calculate the time evolution of the correlation functions 
we have to assume an initial state. The two-photon initial 
state we find under stationary EIT conditions is close to the 
two-polariton ground state. Thus we assume in the following 
as initial state the ground state or a low-temperature thermal 
state. The initial time-independent LL-Hamiltonian for t < 0 
can be diagonalized by means of a Bogoliubov transformation 
which is given by the transformation S = S p { 0) as defined 
above, yielding the Hamiltonian 

H = \p\lllp- (53) 

p^O 


If the initial state is the ground state (T = 0) or a thermal state 
(T > 0 ) of this Hamiltonian the initial correlation functions 
are given by 

= &P,qC0ih(v s \p\/2L T '}, (54) 

where Lt denotes the thermal length corresponding to a finite 
temperature T. 

To begin with, we consider the case T = 0. Using the 
initial correlations ( f54| > in this case and transforming back to 
operators b' p , b p we calculate the time-dependent correlation 
functions, in particular the equal time density-density corre¬ 
lations which are in Luttinger liquid theory universally given 
by equation ( |34] >. The time-dependent correlator G 00 for our 
protocol is now given by 


G^it) = J y (1 - cos (pz)) ((b f p + b- p )(b f _ p + b p )) 

1-00 - a p 

= K{t) / dp-(1 — cos (pz)) 

Jo P 


1 + 


1 - cos(2£(f)/A 0 ) Sin(2£(i)/A 0 ) 


L 2 0 P 2 ~ 1 


V L oP 2 - 1 


(55) 


where we introduced a cutoff a to treat UV divergences. The 
full integral can only be evaluated numerically, but we can 
analyze the limiting cases for small and large momenta p 
giving us asymptotic results for large and small distances z, 
resp. Since A corr (f) is growing only logarithmically with 
A'(f) we assume A corr (f) L 0 . For pL 0 1 the terms 
cos(2£/Ao) an d sin(2£/Ao) oscillate quickly and thus aver¬ 
age to zero. The remaining integral can be calculated analyti¬ 
cally and gives the fully adiabatic correlations 

G = — A'(t) In (^ —+ const. (56) 


Rf’W = ft' 1 " e .sc°/i.) Ri 2) (») (52) 

where £(f) = 2L 1 ^ yJp 2 L\ — 1. It is now straightforward 

to calculate the time evolution of the operators b^ p {t) and b p (t) 
by inverting the transformations. 


leading to a power law in the density-density correlations 
as in the initial ground state but now decaying with the 
much smaller exponent 2K(t) <C 1, thus indicating a quasi¬ 
crystalline order. For pLo -C 1, on the other hand, £(i) be¬ 
comes purely imaginary and approximates to ln(AT(0)/AT(f)) 
such that the last line in ( |55] > approximates to jJjpr cancelling 
















13 



Figure 8. Spatial envelope exp{—2 G^ ! / > (z)} of the CDW correlation 
functions for storage the storage protocol ( |49| ). Comparison of full 
numerical integration of the time evolution using interpolated K(Q) 
to the analytical result {55) using formula {36), where we fix the ini¬ 
tial 0o for both results. The correlation function shows a crossover 
from an adiabatic regime, where the spatial decay fits well to a Gaus¬ 
sian function, to a diabatic regime, showing a power law with ini¬ 
tial exponent 2 Kj. The dashed vertical lines indicate the distance 
2 = L corr . Inset: Plotting the correlations on linear scale show that 
for small distances the correlation functions are well described by a 
Gaussian function with the width a = L corr . 


the prefactor I\ (t) and indicating a power law decay with the 
initial exponent 2/('(()). These limits agree with the crossover 
from an adiabatic to a diabatic regime which we expected 
from the gaplessness of the model. 

We calculate the density-density correlation function 
after the storage numerically and compare it to a full numeri¬ 
cal integration of {45) using an interpolated A'(0) according 
to the DMRG results in Fig. [5] In Fig. [8] we show space- 
dependent amplitude e~ 2G '“ i <;'(Tr of the density-density corre¬ 
lations {34) of both results, for the same intial ©o. 

We plot the distance in units of the length L stor , debited 
by the distance a wave-packet propagates during the storage 
protocol, which we will introduce below. 

The correlations show the expected crossover at length 
scale L C orr from an adiabatic to a diabatic regime, where 
asymptotically the correlation function decays as a power law 
with initial exponent K$. Note that we choose the same ini¬ 
tial ©o for both, analytical and numerical calculation, which 
yields different initial Kq, since the analytic calculations use 
the approximate relation between K and 0. As can be seen 
from the inset in Fig. [8] the adiabatic spatial regime is well 
described by a Gaussian proble exp{— z 2 /2<t 2 }. From bt- 
ting this proble to the calculated correlations we can extract 
the full width at half maximum a of the Gaussian giving the 
value where the amplitude of density-density correlations has 
dropped to 1/2, yielding a length scale for the crossover from 
adiabatic to diabatic regime corresponding to L corr . 


E. Limitations. 

In experimental realizations of the proposed scheme only 
bnite medium lengths are available. Thus the time scale r 
of the protocol has to be limited such that the distance the 
polaritons travel during storage is less than the medium length. 
During the protocol {4 S } the polaritons propagate the distance 

-^stor — j ds v g (s) < Vg(0)@ o K O T. (57) 

As for a bnite system this scale is ultimately limited, we com¬ 
pare L s tor to the correlation length scale L corr . The timescale 
r drops out and we get the best possible L corr for a given L stor 
as 


Lcorr _ 7T PoL a bs |^| j / L .o \ (581 

W “ 45 ~Kq 7 n \K~J ’ ( 

where the initial Kq is also dependent on po°B and 0D|;- A 
feasible Kq can be read off Fig. [5] showing that it is possible 
to get strong correlation functions on the order of the medium 
length for sufficiently small bnal Kf. 


F. Corrections 

In the following we consider corrections to the results due 
to deviations of the initial state from the ground state and ad¬ 
ditional losses during storage due to nonadiabatic couplings. 

Initial excitations. As we have shown numerically the ini¬ 
tial state created under stationary propagation of a two-photon 
pulse into a Rydberg medium is close to the ground-state. As 
the Rydberg polaritons can only be excited inside the EIT 
window CD we can estimate the maximum allowed momen¬ 
tum buctuations by |fc| < |/c max | = f2 2 /|r|c. Excitations 
with a kinetic energy corresponding to k > A: Irlax couple to 
bright-state polariton degrees of freedom which propagate al¬ 
most at the vacuum speed of light, or, in the case of near 
single-photon resonance, are absorbed and thus quickly dis¬ 
appear. Therefore, modeling excitations of the initial state 
within the EIT window by a bnite initial temperature, we 
can estimate an upper bound for this temperature by fee To < 
(fig + g 2 y/n)/2m\Y\c, yielding a thermal length scale HTl of 
(ks = 1 ) 


Lt> ^ = 2p 0 aJ A| 1 


7 tT 


7 ODb 


(59) 


We then take the thermal excitations into account in the ini¬ 
tial correlation functions {54) of the time-dependent Luttinger 
calculation. The thermal length scale Lt marks a crossover of 
the initial state from power law decay of the spatial correlation 
functions to an exponential decay ED. We again calculate 
the time evolution of the density-density correlation functions 
during the storage and plot the resulting functions in Fig. 10 
We bnd that while the correlation function of the thermal ini¬ 
tial state decays exponentially for z L T , the correlation 
function after storage shows adiabatic following for a region 
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Figure 9. Correlation length scales after storage for a thermal initial 
state. The red squares show the thermal length scale from a Gaussian 
fit to the correlation function exp{ — which we compare 

to the geometric mean of Lt and L corr for initial ground state. As 
can be seen, for small thermal length scales the resulting correlation 
length scale is given by this geometric mean and saturates asymptot¬ 
ically tO Lcott for Lt 7^ L C o rr- 



Figure 10. Correlation function exp{— 2G,/,^(z)} for a fixed finite 
initial Temperature To and two different crossover length scales. The 
resulting crossover increases with Lo- Vertical lines correspond to 
the length scales Lq = 21.6 and 


larger Lt- In particular by fitting a Gaussian to the correlation 
function and extracting the crossover length scale, we find that 
after storage instead of the T = 0-correlation length L corr the 
relevant length scale is L CO rr.T which is 


^corr.T — ^ \/ I-corr Lt ■ 


(60) 


Numerically we find that this condition holds as long as Lt < 
2L corr , while for large initial thermal length, L C orr,T asymp¬ 
totically approaches the zero-temperature result L corl . This is 
shown in Fig. [9] 

From this we conclude that if the storage can be done suffi¬ 
ciently slow, if L s tor is sufficiently large, a finite range Wigner 
crystal can be created despite initial thermal excitations or im¬ 
perfections. 

Non-adiabaticity effects. In a cw experiment the validity 
of our model is guaranteed as long as the initial pulse spec¬ 


tral width fits well inside the EIT window defined by con¬ 
dition {nr}. For storage, i.e., a time-dependent control field, 
however, we have to modify condition {TT} and additionally 
consider 6 <C sin0cos0|r e ff |. This restricts the time scales 
on which the input pulse can be stored. For our specific pro¬ 
tocol the nonadiabatic coupling ~ 9 is bounded by its value at 
t = 0. We find 

2K o „ 7 ct 

(K§ - l ) 2 \T\ L abs ' ^ 

Solving for r and using the fact that the maximal allowed r is 
bounded by the condition that the polaritons may maximally 
propagate the system length during storage, i.e., L stoI < L 
we can write this condition as lower bound on the total optical 
depth of the system, 


OD> 


90|F|u g (0) 1 

7t 4 7c 1 — Kq ’ 


(62) 


which obviously diverges for K$ -A 1. However, as we 
can read off from Fig. [5] an initial I\ 0 < 0.8 is experi¬ 
mentally feasible for stationary EIT conditions. Assuming 
|r|u g (0)/7C < 1 we find OD 7> 3, which is easily achiev¬ 
able in current experimental realizations 02- 


SUMMARY 


In summary, we derived an effective model for the interac¬ 
tion of photons in a gas of Rydberg atoms under conditions 
of electromagnetically induced transparency (EIT) in particu¬ 
lar in the regime of large optical depth per blockade distance, 
ODb 7> 1. We showed that under paraxial propagation con¬ 
ditions and for sufficiently low densities of excitations the sys¬ 
tem can be described by an effective model of a single species 
of quasi-particles, called Rydberg polaritons, which becomes 
one-dimensional, if the transverse beam diameter is less than 
the EIT blockade radius. For sufficiently large inter-particle 
distances Rydberg polaritons behave as massive Schrodinger 
particles with repulsive van der Waals type interactions. For 
shorter distances there is a coupling of dark- to fast propa¬ 
gating bright-state polaritons, which gives rise to an effec¬ 
tive loss-mechanism for Rydberg polaritons |[T2 . For an off- 
resonant excitation scheme with finite single-photon detuning 
also bound two-particle states exists. As will be discussed in 
detail elsewhere l27l these states can not be excited for large 
ODb from an initial light field and thus were not considered 
here. We derived conditions where the losses in the effec¬ 
tive Rydberg polariton model are negligible and we can use 
an effective Hamiltonian. The ground-state properties of this 
Hamiltonian were analyzed using numerical DMRG simula¬ 
tions and in terms of a Luttinger liquid model. We showed that 
the regime of strong interactions, quantified by the ratio of in¬ 
teraction to kinetic energy O 7> 1 is very difficult to reach 
under stationary EIT conditions. Increasing the strength of 
the Rydberg interaction leads to an increase of the EIT block¬ 
ade distance, which prevents to reach sufficiently large polari¬ 
ton densities to enter the strong-interaction regime. However, 
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making use of a dynamical slow-down of Rydberg polaritons 
while propagating inside the medium or a storage protocol of 
polaritons in stationary spin excitations allows to decrease the 
kinetic energy by increasing their effective mass without re¬ 
ducing the quasi-particle density. In this way it is possible to 
generate a quasi-crystalline or a charge-density wave (CDW) 
state of stored photons, which is a highly non-classical state 
consisting of an ordered string of single-photon wave-packets. 
This state can be observed by either nonadiabatic release of 


the stored excitations into a train of single photons or by direct 
imaging of the Rydberg ensemble as in ll42l . We analyzed this 
storage in terms of a time-dependent Luttinger liquid model 
and showed that the gapless model leads to a spatial crossover 
in the CDW correlations between an adiabatic regime, exhibit¬ 
ing strong spatial correlations and a diabatic regime, where the 
correlations show the initial power-law decay. 

We thank J. Otterbach, S. Whitlock, M. Weidemiiller and 
H. P. Biichler for valuable discussions. The financial support 
of the DFG through SFB-TR49 is gratefully acknowledged. 
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